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ANALYSIS  AND  SIMULATIONS  OF  OPTICAL  RECTIFICATION  AS  A 
SOURCE  OF  TERAHERTZ  RADIATION 
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I.  INTRODUCTION 

Diode  pumped  laser  systems  based  on  chirped  pulse  amplification  [1]  are  capable  of 
delivering  millijoule  sub-picosecond  pulses  at  repetition  rates  as  high  as  10  kilohertz.  These 
compact  and  efficient  lasers  can  be  used  to  generate  high  average  power  terahertz  (THz) 
radiation  if  an  efficient  photonic  downconversion  scheme  can  be  found. 

Photonic  downconversion  is  a  set  of  methods  whereby  optical  or  infrared  radiation  is 
converted  to  a  lower  frequency  due  to  the  effects  of  the  second  order  susceptibility  present 
in  anisotropic  materials.  With  an  appropriate  pump  source  and  nonlinear  medium  these 
methods  can  be  used  to  generate  THz  radiation  [2,  3,  4,  5,  6,  7,  8,  9,  10,  11].  The  con¬ 
version  efficiency  achieved  by  these  methods  scales  with  several  parameters,  including  the 
pump  intensity,  the  interaction  length,  the  absorption  losses,  the  nonlinearity,  and  the  phase 
mismatch. 

Ultimately,  achieving  breakthrough  conversion  efficiencies  cannot  be  achieved  simply  by 
tuning  parameters  such  as  laser  power  or  interaction  length.  In  fact,  in  the  regime  where 
many  experiments  operate,  the  conversion  efficiency  is  fundamentally  limited  by  the  Manley- 
Rowe  relation  to  a  value  given  by 

Vmr  =  y  (1) 

where  A5  is  the  wavelength  of  the  signal  being  generated  and  Xp  is  the  wavelength  of  the 
pumping  radiation.  In  the  case  of  Ref.  [9],  for  example,  Xp  =  1.06  gm  and  As  =  300  /xm 
so  that  tjmr  =  0.33%.  This  low  value  can  only  be  exceeded  by  designing  a  system  which 
operates  in  a  regime  where  the  assumptions  of  the  Manley- Rowe  relation  do  not  hold.  One 
such  assumption  is  that  the  radiation  consists  of  three  discrete  frequencies.  By  operating  in 
a  regime  where  the  bandwidth  of  the  pump  radiation  is  comparable  to  the  signal  frequency, 
it  may  be  possible  to  obtain  conversion  efficiencies  exceeding  t]mr- 

In  this  report  we  consider  utilizing  ultra-short  laser  pulses  to  efficiently  generate  THz 
radiation  via  a  phase  matched  optical  rectification  process.  Optical  rectification  generates 
radiation  with  a  center  frequency  approximately  equal  to  the  inverse  of  the  pulse  length 
of  the  pump  laser.  If  the  spot  size  of  the  laser  is  many  THz  wavelengths  in  diameter,  the 
process  can  be  phase  matched  by  angle  tuning  in  a  suitable  crystal. 

Manuscript  approved  July  19,  2005. 
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II.  TECHNICAL  BACKGROUND 

A.  Normalizations 

We  use  a  system  of  units  normalized  as  follows.  The  unit  of  time  is  u^1  where  tuT  = 
2tt  x  1012  rad/s.  The  unit  of  length  is  c/cjt,  where  c  is  the  speed  of  light.  The  unit  of  mass 
is  the  electronic  mass,  m,  and  the  unit  of  charge  is  the  electronic  charge,  |e|.  The  unit  of 
density  is  nr  =  muj\/Ane2.  Note  that  this  results  in  the  peculiarity  that  the  unit  of  particle 
number  is  Nt  =  me3/ 4iujTe2.  To  determine  the  value  of  a  normalized  quantity  in  either  SI 
or  gaussian  units,  multiply  the  normalized  quantity  by  the  value  given  in  Table  I.  For  the 
case  of  SI  units,  we  define  the  permittivity  of  free  space,  eo  =  8.85  x  10~12  F/m,  and  the 
impedance  of  free  space,  rj0  =  377  Q. 


B.  Coordinate  System 

Two  natural  coordinate  systems  arise  when  considering  laser  propagation  in  a  crystal. 
In  one  coordinate  system  the  description  of  the  interaction  with  the  crystal  is  simplified, 
while  in  the  other  the  description  of  the  laser  propagation  is  simplified.  The  former  will 
be  described  by  the  “crystal  basis”  defined  by  the  unit  vectors  (u,  v,w).  The  latter  will 
be  described  by  the  “laser  basis”  defined  by  the  unit  vectors  (x,  y,  z).  We  will  consider  the 
laser  basis  to  coincide  with  the  standard  basis. 

In  a  general  coordinate  system,  the  linear  response  of  the  electronic  polarization  P  to  an 
applied  field  E  is  given  by  the  tensor  equation  Pi  =  XijEj  (we  use  the  Einstein  summation 
convention).  By  definition,  Xij  is  diagonal  in  the  crystal  basis.  Furthermore,  for  uniaxial 
crystals  two  of  the  diagonal  elements  are  equal.  By  convention  the  two  equal  elements  are 
Xuu  =  Xw  In  this  case,  radiation  with  a  polarization  component  in  the  w  direction  is  called 
an  extraordinary  wave,  while  radiation  with  no  polarization  components  in  the  w  direction 
is  called  an  ordinary  wave. 

The  laser  basis  is  defined  such  that  the  central  wavevector  of  the  laser  points  in  the  z 
direction.  The  orientations  of  x  and  y  with  respect  to  the  laser  polarization  can  only  be 
made  meaningful  after  specifying  the  relation  between  the  laser  basis  and  the  crystal  basis. 
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TABLE  I:  Normalization 


Quantity 

Symbol 

SI  Unit 

cgs  Unit 

Unit  Value 

Time 

C Orp 

LO'ji 

159  fs 

Length 

cj  COT 

cj  LUX 

47.8  fj,m 

Density 

Tl'p 

e^mu^/e2 

muo\f'pKe2 

1.24  x  1016  cm'3 

Particle  Number 

eo  mc3/ujTe2 

me3  /  4ttu)xc2 

1.35  x  109 

Electric  Field 

Ex 

maux/c 

meux/e 

107  MV/cm 

Magnetic  Field 

mux/e 

mcujx  /  e 

35.7  T 

Current  Density 

rixec 

nxec 

59.5  A/cm2 

Charge  Density 

nxe 

nxe 

0.0020  C/cm3 

Polarization 

tit&c/ujt 

nxec/ujx 

9.5  x  10“8  C/m2 

Susceptibility0 

1 

1/47T 

Susceptibility  (2nd  order) 

1  /  Ex 

1  /  47 r  Ex 

9.3  x  10'11  m/V 

Susceptibility  (3rd  order) 

1  /El 

1/4ttE% 

8.7  x  10'21  m2/V2 

Energy 

9 

me 

o 

mcz 

8.19  x  10'14  J 

Energy  Density 

mc2nx 

mc2nx 

1.02  kJ/cm3 

Power 

mc2cox 

mc2u>x 

0.52  W 

Fluence 

Ex/yo^x 

cEy  / \ixiot 

4.83  J/cm2 

Intensity 

It 

E^/rio 

cE 2  /  4tt 

3.04  x  1013  W/cm2 

Two  Photon  Coefficient 

U>t  /  clx 

lot  /  clx 

0.0069  cm/GW 

“For  the  SI  units,  we  use  the  convention  whereby  €q  is  not  absorbed  into  the  susceptibility;  i.e.,  P  =  cqxE. 


The  laser  basis  and  the  crystal  basis  are  connected  by  a  rotation  operator 


T  =  Ry(-9)Rz(-<j>) 


(2) 


where  Ry  and  Rz  represent  right-handed  rotations  about  the  y  and  z  axes, 
explicitly, 


T  = 


^  cos  9  cos  4>  —  cos  9  sin  <fi  —  sin  9  ^ 
sin  4>  cos  4>  0 

y  sin  9  cos  (j>  —sin  9  sin  $  cos  9  y 


Written  out 


(3) 
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Physically,  the  operator  T  corresponds  to  the  following  procedure.  Position  the  crystal  such 
that  the  crystal  basis  coincides  with  the  laser  basis.  Rotate  the  crystal  —4>  about  the  2:  axis, 
then  rotate  it  —6  about  the  y  axis.  With  this  procedure,  waves  polarized  in  the  ^-direction 
are  extraordinary  waves  while  waves  polarized  in  the  y-direction  are  ordinary  waves. 

Mathematically,  the  operator  T  takes  a  vector  expressed  in  the  crystal  basis  and  gives 
the  same  vector  expressed  in  the  laser  basis: 

A  lb  =  TAqb  (4) 

Here,  the  subscript  CB  refers  to  the  crystal  basis  and  the  subscript  LB  refers  to  the  laser 
basis.  If  the  basis  vectors  themselves  are  expressed  in  the  standard  basis  ,  they  satisfy 

(u,  v,  w)  =  T(x,y,z)  (5) 

An  operator  L  is  transformed  according  to 

Llb  =  TLcbT-'  (6) 

Because  T  is  an  orthogonal  matrix,  its  inverse  can  be  computed  simply  by  taking  the  trans¬ 
pose 

T-i  =  tt  (7) 

C.  Timescale  Separation  and  Laser  Frame 

Numerical  models  of  optical  rectification  can  be  made  more  efficient  by  taking  advantage 
of  the  large  time-scale  separation  between  the  laser  frequency  and  the  THz  frequency.  This  is 
done  by  averaging  over  the  fast  laser  oscillations,  but  not  the  much  slower  THz  oscillations. 
Let  the  real  valued  laser  field  be  denoted  Si  and  the  real  valued  THz  field  by  E{.  The 
complex  amplitude  of  the  laser  field,  <%,  is  then  defined  by 

£i  =  +  c  c  (8) 

2 

The  demand  on  computing  resources  can  be  minimized  by  solving  only  for  the  complex 
envelope  Si  which  varies  much  more  slowly  than  the  real  valued  field  A  similar  notation 
will  be  used  for  all  other  quantities.  For  example,  the  electronic  polarization  is  similarly 
decomposed  into  the  rapidly  varying  part,  denoted  by  the  real  valued  7 \  and  complex  valued 
and  the  slowly  varying  part,  denoted  by  Pi. 
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Further  computational  advantage  can  be  gained  by  carrying  out  the  calculations  in  a 
Galilean  frame  of  reference  moving  at  the  group  velocity  of  the  laser  pulse.  In  this  frame 
of  reference,  the  independent  variables  become  r}  —  z  and  r  =  t  —  ngz.  The  corresponding 
differential  operators  transform  according  to 


9Z  drj  'fl'gdT 

(9) 

dt  =  dT 

(10) 

The  group  index,  often  called  ft,  is  given  by 

dn 

ng  =  n{LJ0)  +  too  ~r 
duj 

U1=LOq 

(11) 

The  utility  of  this  coordinate  system  arises  from  the  fact  that  for  short  pulses  propagating 
in  the  forward  direction  dv  <C  dT.  By  dropping  terms  containing  d ^  from  the  propagation 
equations  one  obtains  equations  that  can  be  integrated  without  the  restrictive  Courant 
condition  that  would  otherwise  apply. 

D.  Boundary  Conditions 

When  a  light  pulse  enters  a  dielectric  the  spatial  length  of  the  pulse  and  the  relative 
strength  of  the  electric  and  magnetic  fields  are  changed.  Suppose  the  dielectric  is  anti¬ 
reflection  coated.  Then  the  energy  in  the  dielectric  can  be  equated  with  the  energy  in  the 
vacuum.  The  energy  of  the  pulse  in  vacuum  is 

Uv  =  J  E2JV  (12) 

where  we  used  Ev  =  BVl  dV  is  a  volume  element,  and  the  integral  is  over  all  space.  The 
energy  in  the  dielectric  is 

u*  =  \j  (El  +  El  +  Ed  ■  P)  dV  (13) 

It  is  easy  to  show  using  Maxwell’s  equations  that  if  the  index  of  refraction  is  n0,  then 

\Bd\  =  n0\Ed\  (14) 

Ud  =  J  n20E2ddV 


Using  P  =  (n2  —  1  )Ed  then  gives 


(15) 
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By  equating  Ud  and  UVi  and  taking  into  account  that  in  the  dielectric  the  pulse  is  spatially 
compressed  by  a  factor  nQ)  we  conclude  that 

El  =  n0Ej  (16) 

B2V  =  %  (17) 

nQ 

It  follows  that  the  Poynting  vector  E  x  B  is  conserved  as  the  pulse  passes  from  vacuum 
into  the  dielectric.  By  expanding  the  Poynting  vector  we  obtain  a  formula  for  the  average 
intensity  in  the  dielectric: 


</>  = 


n0 1  £d  | 


Here  £&  is  the  complex  envelope  described  above.  It  should  be  noted  that  this  intensity  takes 
into  account  the  movement  of  the  energy  associated  with  the  material’s  dipole  moment.  The 
rate  of  power  flow  associated  with  the  fields  only  is  given  by 


(■ffield) 


(i+«r 


E.  Dispersionless  Linear  Propagation  in  a  Uniaxial  Crystal 


Suppose  the  linear  susceptibility  tensor  in  the  crystal  basis  is  given  as 

'  Xo  0  0  ^ 

XCB  =  0  Xo  0 


o  o  Xe 


Then  Eq.  (6)  gives  it  in  the  laser  basis  as 


where 


r  Xu  0  Xc 
Xlb  =  0  Xo  0 

V  Xc  0  X33 


Xn  =  Xo  cos2  9  +  Xe  sin2  0 
Xc  =  (xo  ~  Xe)  cos# sin# 
X33  =  Xe  cos2  6  +  Xo  sin2  9 


7 


Note  that  there  is  no  dependence  on  0,  as  expected  for  a  uniaxial  crystal.  The  exact  wave 
equation  for  the  electric  field  in  a  perfect  insulator  is 

(v2  -  a2)  E  =  «9t2P  +  V  (V  ■  E)  (25) 

Note  that  in  an  anisotropic  medium  the  divergence  term  cannot  be  neglected.  Taking  the 
one  dimensional  limit,  and  using  P*  =  XijEj,  we  obtain  three  equations  in  three  unknowns: 


(a2  -  a2)  ex  =  d\  (XuEx  +  XcEz) 

(26) 

(a2  -  a2)  Ey  =  x0d\Ey 

(27) 

-  d2Ez  =  d2  {xcEx  +  X33 Ez) 

(28) 

The  equation  for  Ez  can  be  solved  immediately  giving 

p  _  p 

1  +  XS3 

(29) 

This  results  in  two  independent  propagation  equations  for  the  transverse  components: 


(d2  -  np2)  Ex  =  0 
(a2  -  n2a2)  Ey  =  0 

where  we  defined 

1  cos2  9  sin2  6 

n2e  n20  +  n2 

nl  =  1  +  Xo 

n2  =  1  +  Xe 

Physically,  n0  is  the  index  of  refraction  for  waves  polarized  in  the  ^-direction  and  uq  is  the 
index  of  refraction  for  waves  polarized  in  the  ^-direction.  Because  na  is  independent  of  0, 
the  y-polarized  waves  are  called  ordinary  waves.  Because  uq  has  an  angular  dependence,  the 
x- polarized  waves  are  called  extraordinary  waves. 

F.  Nonlinear  Polarization  Vector 

The  second  order  nonlinear  polarization  vector  is  related  to  the  electric  field  through  a 
third  rank  tensor  which  is  usually  given  in  the  crystal  basis.  By  writing  out  the  components 
of  the  second  order  polarization  given  by 

pp  +  pW  =  Xijk  (£.  +  £•)  [Ek  +  £k)  (35) 


(30) 

(31) 

(32) 

(33) 

(34) 
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we  obtain  formulas  for  V ^  and  P ^  in  terms  of  the  usual 


( p f  > ) 

^  dn  di2  di3  di4  di5  di6  ^ 

to 

=  T 

C?2i  C?22  (^23  (^24  <^25  G?26 

V  CP  J 

\  ^3i  d3 2  d33  d34  d35  d36  y 

contracted  tensor  elements: 

/  IC.P  ^ 

1C,!2 

Pw|2 

£y£w  T  C.C. 

8y8w  +  c-c- 


(36) 


\  8*8V  +  C.C.  J 


V  ^i2) ) 


f  dn 

di2 

dl3 

di4 

^15 

^16  ^ 

4  T 

d>2i 

^22 

Cl¬ 

io 

CO 

^24 

<^25 

^26 

\  ^31 

^32 

CO 

CO 

"C3 

^34 

<^35 

<^36  j 

8UEU 

8VEV 

£W  Eyj 

f  E  4-  E  E 
8U  Eyj  T  8yj  Eu 


(37) 


y  8UEV  +  8VEU  J 

Given  a  particular  propagation  geometry  the  nonlinear  susceptibility  tensor  can  be  re¬ 
duced  to  a  scalar  denoted  by  deff.  Unlike  the  linear  susceptibility,  des  is  a  function  of  both 
6  and  0  even  for  a  uniaxial  crystal.  The  term  “uniaxial”  is  in  this  sense  a  misnomer  since 
the  crystal  is  only  uniaxial  in  the  linear  approximation.  Calculation  of  de g  by  hand  is  time 
consuming,  but  by  implementing  the  following  procedure  using  symbolic  math  software  it 
can  be  determined  easily.  Consider  first  the  slowly  varying  second  order  polarization  vector. 
Express  the  six  element  column  vector  from  Eq.  (36)  in  terms  of  the  laser  field  components 
expressed  in  the  laser  basis.  For  example,  in  place  of  8U  put 


T~  i 


1  8*\ 

8, 


1 1\ 


(38) 


w  w 

Next,  multiply  out  Eq.  (36)  given  the  contracted  tensor  elements  for  the  crystal  of  interest 
(fortunately  most  of  the  elements  will  usually  vanish).  Also,  it  is  usually  a  good  assump¬ 
tion  that  the  longitudinal  field  components  can  be  neglected  compared  to  the  transverse. 
Next,  determine  which  polarization  components  are  of  interest  given  the  phase  matching 
constraints.  For  type  I  phase  matching  there  are  two  possibilities.  In  the  first  case,  the  THz 
wave  is  an  ordinary  wave  and  Eq.  (36)  will  reduce  to 


d2)  =  deff(M)l4|2 


(39) 
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In  the  second  case  the  THz  wave  is  an  extraordinary  wave  and  Eq.  (36)  will  reduce  to 

pj2>  =  deS(e^)\£y  I2  (40) 

For  type  II  phase  matching,  there  are  two  more  possibilities.  They  are  as  follows: 

pw  =  def!(e,<t>)(£;£v  +  c.c.)  (41) 

Pf>  =deS(9,4>Wx£v  +  c.c.)  (42) 

Note  that  although  we  have  used  the  symbol  de ^  four  times,  in  each  case  the  angular  de¬ 
pendence  is  different.  The  procedure  for  computing  the  value  of  de ff  to  use  for  the  rapidly 
varying  polarization  is  exactly  analogous,  with  the  caveat  that  the  factor  of  4  from  Eq.  (37) 
is  not  absorbed  into  deff. 


III.  DESCRIPTION  OF  SIMULATION  CODE 

A.  Propagation  Equations  for  the  Laser 

For  the  rapidly  varying  time  scale,  we  account  for  dispersion  by  regarding  the  elements 
of  the  susceptibility  tensor  as  operators  of  the  form 


Xij  —  Y2 


(—i)k  dkx 


ij 


k= 0 


du )k 


u  o 


dk 

dt* 


(43) 


which  satisfy  the  equation 

Vi  =  Xij£j  (44) 

Usually  the  operators  are  known  in  the  crystal  basis,  in  which  case  the  laser  basis  operators 
can  be  computed  in  a  way  exactly  analogous  to  the  dispersionless  case: 


Xu  =Xo  cos2  9  +  Xe  sin2  6 

Xc  =  (Xo  ~  Xe)  cos# sin  0 
X33  =  Xe  cos2  6 +  Xo  sin2  9 


(45) 

(46) 

(47) 


The  wave  equation  for  the  complex  laser  field  is  obtained  by  inserting  Eq.  (8)  and  the 
corresponding  form  for  V  into  Eq.  (25).  In  the  one-dimensional  limit  this  gives 


(dz  —  iko)2  —  ( dt  +  ioJo)2  Sx  —  (dt  +  iu>o)2V: 


(48) 
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(dz  -  ik0)2  -  (dt  +  iu0)2  £y  =  ( dt  +  iu o)2Vy  (49) 

Sx  =  -Vz  (50) 

The  polarization  components  are  given  by  the  operator  equations 

Px  =  XuSx  +  Xc£z  +  $x  (51) 

Py  —  Xo£y  +  Sy  (52) 

Pz  —  Xc£x  +  X33  £z  +  (53) 

where  Sz  represents  all  nonlinear  contributions  to  the  polarization. 


Combining  the  equations  for  the  x  and  2  components  involves  an  operator  inversion.  The 
exact  equation  for  the  2-component  is 

(1  +  X.33)£Z  =  ~Xc£x  —  <SZ  (54) 


To  invert  the  operator  on  the  left  hand  side,  first  define 


^33  =  1  +  X3s(w0) 


.  ~  (-j)k  dkX33  dk 

k\  duk  dtk 

K=L  UJ  0 


so  that 


+  <5)  £z  —  ~Xc£x  —  $z 


Far  from  any  resonances,  6  n33.  In  this  case,  iteratively  multiplying  by  the  conjugate  of 

the  operator  on  the  left  gives 

&  =  -C(xA  +  4)  (58) 


where  to  second  order 


c~±- A+it 

^  o  4^6 

n33  n33  n33 


The  propagation  equations  for  the  ordinary  and  extraordinary  waves  are  then 

(dz  —  iko)2  —  ( dt  +  iu o)2  £x  =  (dt  +  iu) o)2  (xn  —  Cxl)£x  —  Cxc^z  +  $x  (60) 
(dz  -  ik0)2  ~  ( dt  +  ^o)2]  £y  =  {dt  +  W o)2(x0£y  +  Sy)  (61) 
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B.  Propagation  Equations  for  the  THz  Wave 


For  the  THz  field,  dispersion  is  not  introduced  by  regarding  the  susceptibility  as  an 
operator,  but  rather  by  coupling  the  field  to  a  population  of  harmonic  oscillators.  The 
propagation  equations  can  be  written  down  immediately  by  taking  the  laser  propagation 
equations,  replacing  the  susceptibility  operators  with  the  corresponding  constants,  taking 
h  =  0,  oj0  —  0,  =  0,  and  by  substituting  Si  +  Hi  for  Si.  The  result  is 

(d2z  -  riffi)  Ex  =  d2  Sz  +  Hz)  +  Sx  +  Hx  (62) 

L  ^33 

(d2z  -  nX)  Ey  =  dl  (S„  +  Hy)  (63) 

where  Si  represents  all  nonlinear  contributions  to  the  polarization,  and  Hi  is  the  polarization 
due  to  the  harmonic  oscillators.  The  harmonic  oscillator  polarization  is  most  conveniently 
computed  in  the  crystal  basis  where  each  component  satisfies  an  independent  equation  of 
the  form 

(dt  +  "A  +  fi,2)  Ht  =  PiEi  (64) 

Here  the  index  i  refers  either  to  the  u,  v,  or  w  coordinate,  pi  is  the  coupling  strength  for 
the  given  polarization,  is  the  damping  constant  for  the  given  polarization,  and  ft,2  is  the 
resonant  frequency  for  the  given  polarization.  Once  Hi  is  found  in  the  crystal  basis  it  can 
easily  be  expressed  in  the  laser  basis  using  H lb  =  THcb- 


Physically,  the  polarization  term  Hi  represents  the  coupling  of  the  THz  radiation  to 
the  vibrations  of  the  crystal  lattice.  In  practical  terms,  however,  what  is  usually  given  is 
the  susceptibility  and  its  derivatives  with  respect  to  frequency.  The  parameters  pi  and 
can  be  computed  from  this  information  by  a  second  order  matching  of  dispersion  relations. 
Suppose  the  dispersion  information  for  ordinary  waves  is  given  by  a  susceptibility  x  and  its 


(65) 

(66) 

(67) 


Here,  the  prime  denotes  differentiation  with  respect  to  angular  frequency  with  evaluation  at 
cjs.  The  parameter  cus  can  be  any  frequency  in  the  THz  range  where  x  and  its  derivatives 
are  known.  The  case  for  Xe,  and  pw  is  exactly  analagous. 
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C.  Numerical  Solution  of  the  Propagation  Equations 

1.  Time  Centering  and  Discrete  Grid 

Numerical  solution  of  the  propagation  equations  is  carried  out  by  differencing  the  propa¬ 
gation  equations  on  a  discrete  grid.  The  “spatial”  grid  corresponds  to  the  coordinate  r  and 
the  “time”  levels  correspond  to  the  coordinate  77.  The  electric  field  vectors  Ei  and  Si  are 
evaluated  at  the  same  spatial  grid  points  (i.e.,  a  staggered  spatial  grid  is  not  used).  At  the 
beginning  of  a  simulation  cycle,  the  laser  fields  Si  are  known  at  time  levels  n  and  n  +  1.  The 
THz  fields  Ei  are  known  at  time  levels  n  +  1/2  and  n  +  3/2.  The  simulation  cycle  consists 
of  advancing  the  laser  fields  to  time  level  n  +  2  and  the  THz  fields  to  time  level  n  -f  5/2. 
The  laser  fields  are  advanced  first.  The  laser  fields  from  time  level  n  are  used  to  implicitly 
advance  the  linear  part  of  the  propagation  equation.  The  laser  fields  from  time  level  n  +  1 
and  the  THz  fields  from  time  levels  n  +  1/2  and  n  +  3/2  are  used  to  include  the  nonlinear 
part  as  an  explicitly  known,  centered,  source  term.  The  THz  fields  are  advanced  second. 
The  THz  fields  at  n  +  3/2  and  the  laser  fields  at  n  +  2  are  used  in  the  explicit  advance  of 
the  THz  fields  to  time  level  n  +  5/2. 


2.  Laser  Propagation 


Expanding  the  equations  for  the  laser  evolution  results  in  a  large  number  of  terms.  To 
handle  this  complexity  we  utilize  symbolic  math  software  to  automatically  generate  the  code 
needed  to  advance  the  Eqs.  (60)  and  (61).  The  procedure  utilized  by  the  software  program 
is  as  follows.  First,  make  the  substitutions  dz  —  dv  —  ngdT ,  dt  —  dT ,  and  ko  =  cj0?V  Next, 
expand  the  equations  throwing  out  all  derivatives  higher  than  the  first  order  in  77  and  the 
second  order  in  r.  Next,  substitute  difference  operators  for  the  differential  operators  as 
follows: 


dX 


[S(i  +  1,  n)  +  S(i  +  1,  n  +  2)J  -  [S(i  -  1,  n)  +  S{i  -  1,  n  +  2)] 

4Ar 


(68) 


[S{i  +  1,  n)  +  S (i  +  1,  n  +  2)]  —  2 [£(i,  n )  +  £{i ,  n  +  2)]  +  [S(i  —  1,  n)  +  £{i 

2A  t2 


a  c  ^  [£{i  +  l,n  +  2)  -  S(i 

OnrL  * 


1,  n  +  2)]  —  [S(i  +  1,  n)  —  S(i 
4Ar)Ar 


!.«)] 


1,  Ti  +  2)] 

(69) 

(70) 
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dv£ 


£(i,n  +  2)  -  £(z,n) 
2Ar? 


(71) 


Here  £  is  either  electric  field  component,  z  is  the  grid  cell,  n  is  the  “time”  level,  At  is  the 
cell  size  and  A rj  is  the  timestep.  The  resulting  equation  will  always  take  the  form 


T\£[i  —  1, 7i  T  2)  +T2£  (z,  77, -j- 2)  -\-T$£  (z  4- 1, 77.4-2)  —  K\£(i  1,  tl)  4~  K2£  (z,  n)  4~  K^£  (z  T 1,  n)  +  C7 

(72) 

Note  that  the  nonlinear  terms  are  contained  in  C  and  are  centered  since  they  are  known 
at  time  level  n  -f  1.  The  final  outputs  from  the  symbolic  math  software  are  expressions  for 
Tl,  Ku  and  C.  To  eliminate  transcription  error,  it  is  best  to  have  the  software  generate  the 
expressions  in  a  format  compatible  with  the  programming  language  being  used  for  the  sim¬ 
ulation  so  that  the  expressions  can  be  pasted  directly  into  the  code.  Once  these  expressions 
are  in  the  code  a  standard  tri diagonal  algorithm  can  be  used  to  advance  £{. 


3.  THz  Propagation 


Rewriting  the  THz  propagation  equations  in  terms  of  the  laser  frame  coordinates  and 
taking  — >  0  gives 


(vedv  +  dT)  Ex  = 


dT 


ng  ~  ne 


-rr{Sz  +  Ht)  +  sx  +  Hx 


n 


(■ v0dv  +  dT)  Ey  = 


33 


dT 


n„  —  ni 


[S'y  4-  Hy 


where 


2  n, 


ve  = 


ne  -  nl 


2  n. 


vn  = 


A  ~  nl 

The  equation  for  the  lattice  vibration  in  the  laser  frame  is 


(73) 

(74) 

(75) 

(76) 


(aT2  +  VidT  +  n2)  ^  =  PiEi  (77) 

The  nonlinear  polarization  Si  is  computed  from  the  laser  field  using  Eq.  (36). 

The  problem  is  to  advance  Ej  from  timestep  n  +  3/2  to  timestep  n  +  5/2  given  that  we 
know  Hi  at  timestep  n  +  3/2  and  at  timestep  n  +  2.  Since  Hi  is  effectively  a  function  of 
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Ei ,  the  equations  for  Ei  can  be  put  in  the  form  of  the  flux  conservative  initial  value  problem 
d^Ei  =  — drFi(Ei )  where 


Ex  Sx  +  Hx  —  Xc(Sz  +  Hz)/n\z 
vo  ve(n2g  -  n2e) 


(78) 


_  ^v_  _  +  Hy 

V  v0  v0(n2g  -  n20) 

These  can  be  solved  using  the  standard  Lax-Wendroff  method  where  Hi  is  regarded  as  a 
function  of  Ei  and  Si  is  regarded  as  constant. 


4-  Longitudinal  Fields 

Although  the  longitudinal  fields  do  not  appear  explicitly  in  the  propagation  equations, 
they  are  needed  to  compute  the  nonlinear  source  terms  and  the  lattice  vibrations.  In  the 
laser  advance,  the  nonlinear  source  terms  Si  are  needed  at  time  level  n  +  1,  and  therefore 
all  the  components  Ex  and  Ex  are  also  needed  at  time  level  n  - b  1.  These  will  be  available 
provided  all  the  components  including  £z  were  advanced  to  time  level  n  +  2  on  the  previous 
cycle.  To  accomplish  this,  we  first  advance  8X  as  described  above,  and  then  use  the  linearized 
Eq.  (58)  to  estimate  £z : 

£z  «  -Cxctx  (80) 

In  the  THz  advance,  the  lattice  vibrations  Hx  depend  on  the  longitudinal  field  Ez.  Therefore, 
at  the  midpoint  of  the  Lax-Wendroff  calculation  we  estimate  Ez  based  on  the  slowly  varying 
analog  of  Eq.  (58): 

r- ,  _  XcEx  +  Hz  T  Sz 
E z  1 

1  +  X33 

We  then  recalculate  Hx  based  on  this  new  field.  At  the  end  of  the  Lax-Wendroff  calculation, 
we  repeat  this  procedure  twice  in  order  to  iteratively  improve  the  accuracy  of  both  Ez  and 
Hx  for  use  during  the  next  cycle. 


5.  Magnetic  Field  and  Energy  Diagnostic 

Although  the  magnetic  field  is  not  needed  to  advance  the  propagation  equations,  it  is 
useful  for  diagnostic  purposes.  In  particular,  it  can  be  used  in  the  calculation  of  the  energy 
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contained  in  the  laser  oscillations  and  the  THz  oscillations  so  that  the  conversion  efficiency 
can  be  calculated,  and  energy  conservation  can  be  verified. 

For  the  magnetic  field  associated  with  the  laser  energy,  we  first  compute  the  linear  part 
of  the  polarization  using  Eq.  (44)  and  then  use 


Bx=  £y  +  Vy 

(82) 

Tig 

£r  H-  Vr 

By=  X  X 

(83) 

ng 

Bz  =  0 

(84) 

which  is  an  approximation  based  on  Ampere’s  law.  Note  that  dispersion  can  affect  the 
proportion  of  energy  carried  by  the  magnetic  field.  Remembering  the  summation  convention, 
the  energy  density  carried  by  the  fast  laser  oscillations  is 


uf 


£j£*  +  BjB*  t  EtVi  +  EjVt 


(85) 


For  the  magnetic  field  associated  with  the  THz  energy,  the  linear  part  of  the  polarization 
can  be  found  from 

Px  =  (nj  -  1  )EX  +  Hx  (86) 

Py  =  ( no  ~  l)  By  +  Hy  (87) 

Pz  =  -Ez  (88) 


The  THz  magnetic  field  can  then  be  computed  from  expressions  exactly  analogous  to  those 
for  the  laser  field.  The  energy  density  can  then  be  computed  using 


EiEi  +  BiBi  +  EiPi 


(89) 


IV.  ANALYSIS  OF  OPTICAL  RECTIFICATION 
A.  Driving  Terms  for  Gallium  Selenide 

To  date,  the  conversion  of  optical  to  THz  power  has  been  most  efficient  when  Gallium 
Selenide  (GaSe)  was  used  as  the  nonlinear  material  [9].  We  therefore  use  GaSe  as  an  example 
in  all  that  follows.  The  symmetries  of  the  GaSe  lattice  require  that  all  components  of  the 
nonlinear  susceptibility  vanish  except  for  die,  d2i ,  and  d2 2-  Furthermore,  d16  =  d2i  =  -d22  « 
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50  x  1CT12  m/V.  Using  the  technique  outlined  in  section  II F,  the  transverse  components 
of  the  slowly  varying  nonlinear  polarization  vector  can  be  determined  to  have  the  following 
dependence  on  the  applied  laser  field  and  the  crystal  orientation: 

=  diQCOsd  —  cos2  #sin30|£r|2  +  sin30|£y|2  +  cos9  cos3cf)(£x£*  +  £x£y)  (90) 

Py2^  =  di6  cos2  9 cos 3</)\£x | 2  —  cos30|£0|2  +  cos  0  sin  30(£x£*  -fi  £x£y)  (91) 

Typically,  because  of  phase  matching  considerations,  only  one  of  the  six  driving  terms  from 
the  above  two  equations  will  be  important.  Furthermore,  it  is  virtually  always  the  case  that 
in  order  for  phase  matching  to  be  achieved,  the  driving  term  must  contain  the  polariza¬ 
tion  component  orthogonal  to  the  driven  term.  Hence,  the  expressions  for  the  polarization 
components  can  usually  be  reduced  to  one  of  the  following: 


Px2^1  =  die  cos6*sin30|£y|2 

(92) 

P^'1  —  di6  cos2  0cos30|£x|2 

(93) 

Pj2),n  —  die  cos2  9  cos  30(£x£*  +  £x£y) 

(94) 

Py(2),n  =  cos#sin30(£x£*  +  £*£y ) 

(95) 

Here,  the  superscripts  I  and  II  refer  to  type  I  and  type  II  phase  matching,  respectively.  The 
definition  of  type  I  phase  matching  is  that  the  driving  term  contains  only  one  polarization 
component.  The  definition  of  type  II  phase  matching  is  that  the  driving  term  contains  both 
polarization  components. 

In  the  case  of  the  nonlinear  polarization  at  the  laser  frequency,  the  expressions  for  the 


transverse  components  are 

Vx2^  =  4di6  cos2  9  cos  30  ( £xEy  +  £yEx )  —  cos  9  sin  30  ( £XEX  cos2  9  —  £yEy^  (96) 

Py2^  =  4 d16  cos  30  (£XEX  cos2  9  -  £yEy^  +  cos  9  sin  30  {£xEy  +  £yEx)  (97) 

For  type  I  phase  matching,  these  can  be  reduced  to  one  of  the  following: 

V^'1  =  4di6  cos2  9  cos  3 <p£xEy  (98) 

V^'1  =  4di6  cos  0  sin  30£yi£E  (99) 
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B.  Phase  Matching  Condition  and  THz  Evolution 


An  expression  for  the  signal  produced  by  optical  rectification  in  the  presence  of  an  ar¬ 
bitrary  phase  mismatch  can  be  found  if  pump  depletion  and  dispersion  are  neglected.  In 
the  case  of  a  type  I  interaction,  dropping  the  dispersive  terms  from  Eq.  (73),  and  assuming 
Sz  «  Sx,  gives 

M„  +  dr)  Ex  =  -pz__dr\£y\2  (100) 

ng  ~  ne 

where  for  a  GaSe  crystal 

<2eff  =  die  cos  9  sin  30  (101) 

Making  a  second  coordinate  transformation  Z  =  r)  and  T  =  t  —  rj/ve  gives 


dzEx  = 


d, 


eff 


2  n. 


dT\Sy 


The  integral  representation  of  Ex  is 


(102) 


The  integral  cannot  be  solved  in  this  form  because  £y  is  a  function  of  both  T  and  Z.  However, 
if  pump  depletion  and  dispersion  are  neglected,  Ey  is  independent  of  77.  We  therefore  rewrite 
the  integral  as 

Ex  =  (104) 

which  can  be  solved  immediately.  Defining  e  =  rj/ve  we  obtain 

Ex  =  7  (IWI2  -  l^T  -  e)i2)  (105) 

For  finite  e,  as  the  propagation  distance  77  increases  the  two  terms  in  Eq.  (105)  become  well 
separated  in  time  and  the  peak  signal  amplitude  converges  to  the  finite  value 

1 T  _  ^eff|^o|2  ^ 

-Epeak  —  ,  2  or  V^>ve^L  (106) 

\ne  ~  ng\ 

where  £0  is  the  peak  electric  field  of  the  laser  pulse  and  rL  is  the  laser  pulse  length.  Con¬ 
vergence  of  the  signal  to  a  finite  value  is  not  characteristic  of  a  phase  matched  interaction. 
The  case  where  e  =  0  can  be  analyzed  by  expanding  Eq.  (105)  in  powers  of  e: 


1 -1  deft 

bx  = - 77 

2nn  ' 


tdl  e2dl 


dx  +  ir  + 


3! 


+  . . 


\£yij)\2 


(107) 
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Signal  Wavelength  (p m) 


FIG.  1:  Angle  tuning  curve  for  optical  rectification  in  GaSe  using  800  nm  radiation  as  the  pump. 

Substituting  e  =  0  then  gives  a  signal  amplitude  that  grows  indefinitely  with  rj.  The  phase 
matching  condition  is  therefore  e  —  0.  The  waveform  of  the  signal  in  this  case  is  just  the 
derivative  of  the  laser  envelope.  For  a  typical  laser  pulse  shape,  the  signal  will  be  a  single 
cycle  pulse  with  center  frequency  us  =  7 t/tl- 

A  more  physically  revealing  form  of  the  phase  matching  condition  is 

ne(us)  =  ng(ujo)  (108) 

In  other  words,  the  phase  velocity  of  the  THz  radiation  must  equal  the  group  velocity  of 
the  laser  pulse.  This  condition  arises  because  the  phase  of  the  source  (i.e.,  the  nonlinear 
polarization  wave)  at  a  given  point  is  determined  by  the  derivative  of  the  envelope  of  the 
laser  pulse  at  that  point.  Therefore  the  phase  of  the  source  stays  synchronous  with  the 
phase  of  the  THz  wave  when  the  laser  envelope  moves  at  the  THz  phase  velocity.  Fig.  1 
shows  the  phase  matching  angle  in  GaSe  as  a  function  of  signal  wavelength  for  an  800  nm 
pump  laser.  The  curve  was  calculated  by  solving  ng(cus)  =  ng(u> 0)  using  the  GaSe  dispersion 
relation  from  Ref.  [12]. 

For  a  perfectly  phase  matched  interaction,  the  signal  amplitude  can  always  be  increased 
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FIG.  2:  Time  dependence  of  the  THz  electric  field  after  1  cm  of  propagation  for  two  matching 
conditions.  The  solid  line  is  the  THz  field  for  a  negligible  phase  mismatch  of  0.01%.  The  dashed 
line  is  the  THz  field  for  a  mismatch  of  2.0%.  The  dotted  line  shows  the  magnitude  of  the  complex 
envelope  of  the  laser  field. 

by  increasing  the  interaction  length,  r).  In  the  case  of  a  significant  phase  mismatch,  Eq.  (106) 
suggests  that  Epeak  could  still  become  large  if  either  |£o|2  or  de q  are  large.  Unfortunately,  as 
discussed  in  section  IV  D,  |£0|2  is  limited  by  the  two-photon  absorption  process.  On  the  other 
hand,  media  with  very  large  values  of  des,  such  as  electric  field  biased  AlGaAs  quantum  wells 
and  GaAs/AlGaAs  asymmetric  quantum  wells  [13,  14],  might  make  it  possible  to  generate 
large  amplitude  THz  waves  without  the  need  for  phase  matching. 

As  an  example  of  a  phase  matched  interaction,  consider  the  generation  of  1.0  THz  radi¬ 
ation  in  GaSe  using  a  pump  intensity  of  1  GW/cm2.  The  signal  frequency  fixes  the  laser 
pulse  length  to  tl  =  500  fs.  Figure  2  shows  the  time  dependence  of  the  THz  electric  field  as 
computed  from  Eq.  (105)  evaluated  after  1  cm  of  propagation  for  two  matching  conditions. 
The  solid  line  corresponds  to  0  =  28.1°  which  gives  |n5(w0)  -  no(us)\/ng(uo)  =  0.01%.  The 
dashed  line  corresponds  to  6  =  20.0°  which  gives  |n5(tu0)  -  ne{uJs)\/ng(LUo}  =  2.0%.  The 
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FIG.  3:  Instantaneous  frequency  shift  Acj/cjo  of  the  pump  after  1  cm  of  propagation  in  GaSe  in 
the  case  where  dispersion  is  neglected.  The  intensity  of  the  pump  is  1  GW/cm2. 

dotted  line  shows  the  laser  envelope  on  a  different  scale.  As  expected  from  the  expansion 
solution,  the  matched  case  is  approximately  the  time  derivative  of  the  Laser  envelope.  In  the 
mismatched  case,  the  time  separation  between  the  positive  and  negative  polarity  portions  of 
the  pulse  increases  with  propagation  distance,  but  the  signal  amplitude  remains  constant. 


C.  Pump  Depletion  and  Pulse  Compression 


The  effect  of  the  signal  field  on  the  laser  is  to  deplete  the  laser  energy  and  compress 
the  laser  pulse.  Assuming  perfect  type  I  phase  matching  (n0  =  uq  =  ng)  and  neglecting 
dispersion  allows  the  propagation  equation  for  the  laser  to  be  reduced  to 


(j-^ri  2'figdyj-  2iu>0ngdvJ  Sy  (^r  ~b  ivo)  <5 y 
Assuming  dT  cj0,  XcSz/n 33  <SX,  and  d2  — >  0,  this  can  be  further  simplified  to 


(109) 
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A  first  order  perturbation  analysis  of  this  equation  gives 

£y  =  £o-—r)cU£0Ex  (111) 

ng 

where  £0  is  the  initial  laser  field.  Inserting  Eq.  (107)  for  Ex  and  taking  e  =  0  gives 


Sy  =  E o  1  +  -4dl^2dT\EQ 


nt 


(112) 


The  argument  of  the  factor  in  parenthesis  gives  the  instantaneous  phase  shift.  Taking  the 
time  derivative  of  this  phase  gives  the  instantaneous  frequency  shift 

Au 


—  =  v2dhd2T\s0 


UJq 


nt 


(113) 


The  instantaneous  frequency  shift  is  plotted  in  Fig.  3  for  the  same  parameters  that  were  used 
in  Fig.  2.  The  frequency  shift  is  proportional  to  the  second  derivative  of  the  laser  envelope. 

It  is  instructive  to  consider  the  frequency  shift  just  derived  in  terms  of  an  analogy  with 
self  phase  modulation.  In  self  phase  modulation,  the  nonlinear  polarization  is  P ^  =  C|£|2£, 
where  £  is  any  component  of  the  laser  field  and  C  is  a  factor  independent  of  £.  This  leads  to 
an  intensity  dependent  refractive  index  n3  =  ni+(7|£|2/2n1,  where  n\  is  the  linear  part  of  the 
refractive  index.  This  leads  to  a  frequency  shift  proportional  to  dT\£\2 .  In  the  present  case, 
the  nonlinear  polarization  is  P O  =  4 de^E£  =  K(dT\£\2)£ ,  where  K  is  a  factor  independent 
of  £.  The  intensity  dependent  refractive  index  in  this  case  is  77-2  =  n\  +  K(dT\£\2)/2n\. 
Comparing  ri2  and  n3  makes  it  obvious  that  if  n3  leads  to  a  frequency  shift  proportional  to 
dr\£\ 2,  then  n2  leads  to  a  frequency  shift  proportional  to  d2\£\2. 

A  photon  kinetic  picture  [15,  16]  can  be  used  to  show  that  in  the  presence  of  group 
velocity  dispersion,  the  frequency  shift  induced  by  the  THz  field  causes  a  portion  of  the 
laser  pulse  to  compress.  Let  r'  =  t  —  ngz'  where  z'  is  the  coordinate  of  a  photon  in  the 
Lagrangian  (as  opposed  to  Eulerian)  picture.  The  motion  of  the  photon  is  then  described 
by 

dP 

-~=(52Au  (114) 

dr] 

where  H( ca0  +  Aw)  is  the  energy  of  the  photon  and  f32  is  the  usual  group  velocity  dispersion 
coefficient  defined  as  the  second  order  term  in  the  sequence 
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FIG.  4:  Length  of  GaSe  crystal  needed  to  maximally  compress  an  800  nm  wavelength  pump  pulse. 


After  the  pulse  propagates  a  distance  L,  a  photon  initially  at  r'  =  0  moves  to 


J  =  [  /32Au>dr] 
Jo 


(116) 


We  approximate  Acj  by  evaluating  Eq.  (113)  at  r  =  0  and  using  for  the  initial  pulse  shape 


S0  =  Sqo  exp  —2  In  2 


(117) 


This  procedure  assumes  that  r'  <C  Tz,.  In  this  case,  the  photon  position  is 

8  In  2  p2uQ  d2eff£$0  r3 

"'TIT' 


(118) 


The  sign  of  (52  is  usually  positive,  which  means  that  the  downshifted  photons  at  the  center 
of  the  pulse  will  tend  to  drift  toward  the  upshifted  photons  at  the  front  of  the  pulse.  An 
estimate  of  the  point  of  maximum  compression  can  be  obtained  by  setting  r'  =  — r^/2  and 
solving  for  L.  This  gives 
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The  compression  length  as  a  function  of  pump  intensity  is  shown  in  Fig.  4  for  the  same 
parameters  that  were  used  in  Fig.  2.  For  a  GaSe  crystal,  (32  ~  4  x  10-16  seconds  for  800  nm 
radiation. 


D.  Two  Photon  Ionization 


The  highest  intensity  useful  for  optical  rectification  is  reached  when  the  electron  hole 
plasma  becomes  dense  enough  to  reflect  THz  radiation.  It  can  be  shown  that  the  absorption 
of  radiation  due  to  ionization  processes  and  carrier  scattering  is  described  by 

f)  T 

T  =  -v*l£-nefThI  (120) 

where  Ui  is  the  energy  required  to  bring  an  electron  from  the  valence  to  the  conduction 
band,  ne  is  the  free  electron  density,  and  a h  is  the  cross  section  for  hole  scattering.  Hole 
scattering  leads  to  collisional  losses  through  inverse  bremstrahlung.  For  GaSe,  the  bandgap 
energy  is  Eg  =  2.0  eV,  and  w  5  x  10-17  cm2  [17].  For  800  nm  light,  the  photon  energy  is 
about  1.55  eV,  so  ionization  occurs  via  a  two  photon  mechanism,  and  Ui  =  3.1  eV.  In  this 
case,  the  free  electron  density  evolves  according  to 


dne  (3  r2 


dr 


(121) 


where  (3  is  the  usual  two-photon  absorption  coefficient.  For  GaSe,  (3  w  0.558  cm/GW  for  800 
nm  radiation  and  a  1  ps  pulse  length  [18].  The  maximum  intensity  useful  for  the  generation 
of  THz  radiation  is  the  intensity  for  which  the  free  electron  density  reaches  the  critical 
density 

■2  1  *  '  ~  '  "  1  (122) 
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where  uj  is  the  frequency  of  the  THz  radiation,  m *  is  the  effective  mass  of  an  electron,  m*h  is 
the  effective  mass  of  heavy  holes,  m*  is  the  effective  mass  of  light  holes,  and  a  is  the  fraction 
of  holes  that  are  heavy  holes.  For  GaSe,  m*  ~  0.2.  To  obtain  an  estimate  of  the  maximum 
useful  intensity,  we  combine  Eqs.  (121)  and  (122),  take  the  laser  pulse  length  to  be  500  fs, 
assume  awl  and  m*h  >  m*: 

'  u>2m*Ui\112 


Imax 


prL  J 


=  2.1  GW/crn 


(123) 
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TABLE  II:  Simulation  Parameters 


Quantity 

Value 

Crystal 

GaSe 

Dispersion  Parameters 

See  Ref.  [12] 

Lattice  damping  constant 

v  —  125  GHz 

Two  photon  ionization  coefficient 

0  =  0 

Laser  Intensity 

1  GW/cm2 

Laser  Wavelength 

Ao  =  800  nm 

Laser  Pulselength 

tl  =  500  fsec 

Signal  Frequency 

1.0  THz 

Phase  Matching  Configuration 

Type  I 

Phase  Matching  Angle 

9  =  28.1° 

Azimuthal  rotation 

-e- 

II 

CO 

o 

o 

Laser  Polarization 

ordinary 

THz  Polarization 

extraordinary 

It  should  be  noted  that  at  this  intensity  the  electron-hole  plasma  will  significantly  modify 
the  dispersion  relation  for  the  THz  radiation  which  will  in  turn  modify  the  phase  matching 
condition.  In  addition,  the  value  of  7max  given  in  Eq.  (123)  will  be  substantially  reduced  if 
the  value  of  /3  measured  in  earlier  experiments  [19,  20]  is  used. 

V.  SIMULATIONS 

The  simulation  model  described  in  section  III  extends  the  regime  of  validity  of  the  anal¬ 
ysis  in  that  it  self-consistently  includes  dispersion  of  both  the  THz  and  laser  pulses,  pump 
depletion,  and  THz  losses  due  to  vibrational  damping.  In  this  section,  we  compare  the  anal¬ 
ysis  with  simulation  results  obtained  using  different  sets  of  approximations.  The  parameters 
used  are  given  in  Table  II. 

First,  the  analytical  result  shown  in  Fig.  2  is  reproduced  by  running  the  simulation  with 
dispersion  neglected.  The  simulation  results  are  shown  in  Fig.  5.  The  two  results  are  almost 
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FIG.  5:  Simulation  with  dispersion  turned  off  showing  agreement  with  the  analytically  predicted 
THz  waveforms  in  the  case  of  a  phase-matched  and  non-phase-matched  interaction  (cf.  Fig.  2). 

indistinguishable.  Note  that  pump  depletion  was  not  neglected,  but  its  effect  was  minimal 
due  to  the  low  conversion  efficiency  over  1  cm  of  propagation. 

From  the  same  simulation,  a  comparison  can  be  made  with  the  analytical  prediction  for 
the  frequency  shift  from  Fig.  3.  The  corresponding  simulation  result  is  shown  in  Fig.  6.  The 
simulation  data  is  displayed  as  a  Wigner  transform  of  the  laser  electric  field.  The  Wigner 
transform  is  characterized  by  the  following  properties:  (i)  when  integrated  over  time  it  gives 
the  frequency  spectrum  of  the  pulse  (ii)  when  integrated  over  frequency  it  gives  the  time 
domain  representation  of  the  pulse  (iii)  it  becomes  a  photon  distribution  function  in  the 
limit  of  geometric  optics.  The  frequency  shift  indicated  by  the  Wigner  transform  is  similar 
to  the  instantaneous  frequency  shift  from  Fig.  3. 

As  an  additional  check,  Fig.  7  shows  the  energy  gained  by  the  THz  radiation  and  the 
energy  lost  by  the  laser  pulse  for  the  same  simulation  used  to  generate  Figs.  5  and  6.  The 
simulation  conserves  energy  to  a  high  degree  of  accuracy. 

Next,  the  simulation  is  repeated  with  dispersion  included,  and  the  propagation  distance 
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FIG.  6:  Wigner  transform  of  the  laser  electric  field  after  1  cm  of  propagation  in  the  absence  of 
dispersion. 


FIG.  7:  Energy  conservation  in  the  case  of  zero  dispersion. 

is  extended  far  enough  to  observe  saturation  of  the  conversion  efficiency.  The  damping 
constant  associated  with  the  lattice  vibrations  is  set  to  v  =  125  GHz  in  order  to  give  an 
absorption  coefficient  consistent  with  Ref.  [6].  The  conversion  efficiency  as  a  function  of 
propagation  distance  is  shown  in  Fig.  8.  A  conversion  efficiency  of  1.5%  is  obtained  before 
saturation  occurs  at  about  z  —  5  cm.  Note  that  the  Manley-Rowe  limit  is  exceeded  by  about 
a  factor  of  5.  The  THz  waveform  at  z  =  5  cm  is  shown  in  Fig.  9.  The  waveform  departs 
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FIG.  8:  Energy  conservation  with  dispersion  and  vibrational  damping. 

from  the  single  cycle  pulse  characteristic  of  optical  rectification  because  of  the  dispersive 
terms  Hi. 

Finally,  Fig.  10  shows  the  laser  pulse  shape  at  four  propagation  distances.  As  localized 
frequency  shifts  develop,  pulse  compression  occurs  at  the  head  of  the  pulse  due  to  group 
velocity  dispersion.  In  this  case,  the  point  of  maximum  compression  occurs  at  about  z  =  2.1 
cm,  which  is  about  twice  the  analytical  prediction  from  Fig.  4.  The  fact  that  the  group 
velocity  and  pulse  shape  both  change  during  the  course  of  the  interaction  causes  the  phase 
mismatch  to  increase  with  increasing  z.  This  eventually  leads  to  saturation  of  the  conversion 
efficiency.  It  should  be  noted  that  in  this  example  the  two-photon  ionization  rate  increases 
by  about  a  factor  of  12  as  the  laser  propagates  from  z  =  0  cm  to  z  —  2  cm.  At  the  same 
time,  the  pulse  length  is  is  reduced  by  a  factor  of  5.  The  electron  hole  plasma,  therefore,  is 
about  twice  as  dense  at  2  =  2  cm  than  at  z  —  0  cm.  This  could  strongly  affect  the  results, 
and  will  be  the  subject  of  future  study. 

VI.  SUMMARY 

Optical  rectification  of  short  laser  pulses  in  nonlinear  crystals  can  generate  high  peak 
power  THz  radiation.  The  GaSe  crystal  appears  to  be  promising  for  this  application  due  to 
its  high  nonlinearity,  low  losses,  and  favorable  phase  matching  characteristics.  Analysis  and 
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FIG.  9:  Time  domain  waveform  of  the  THz  signal  after  5  cm  of  propagation 

simulations  both  indicate  that  an  800  nm  laser  pulse  with  a  pulse  width  of  500  fs  can  be 
phase  matched  in  GaSe  to  produce  1.0  THz  radiation  with  a  peak  power  of  about  1  MW.  The 
conversion  efficiency  saturates  at  about  1  percent  due  to  frequency  shifts  in  the  laser  pulse 
which  cause  it  to  compress  and  accelerate.  Also,  it  should  be  noted  that  the  configuration 
considered  here  requires  that  index  matching  materials  be  used  to  couple  radiation  into 
and  out  of  the  GaSe  crystal.  In  the  future,  the  electron-hole  plasma  will  be  incorporated 
directly  into  the  simulation  model  so  that  its  effects  on  the  conversion  efficiency  can  be  more 
accurately  determined. 
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FIG.  10:  Temporal  shape  of  the  laser  pulse  evaluated  at  four  propagation  distances.  Note  that  the 
increasing  laser  intensity  will  increase  the  two- photon  ionization  rate. 
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